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Preface 


This  work  is  a  continuation  of  work  done  by  Shell  in  the  stability  evaluation 
of  a  scalar  control  applied  to  a  time  periodic  system.  The  scalar  controlled  satellite 
is  examined  in  a  new  controlled  modal  coordinate  system.  This  coordinate  system 
decouples  the  four  modes  of  the  linearized  system  with  the  controller.  The  stability 
of  each  mode  can  easily  be  evaluated  in  the  uncoupled  form.^  This  technique  is  . 

■  ■‘V 

demonstrated  on  a  well  defined  problem  and  can  easily  be  extended  to  higher  order  _  ^ 

systems.  The  comparison  of  the  linear  and  non-linear  equations  in  the  controlled  f  ‘  1 

modal  coordinates  demonstrates  some  of  the  limitations  of  linearized  equations  in 
this  type  of  non-linear  system. 

I  must  thank  Dr.  Wiesel  for  his  consultation  on  the  computer  programs  and 
help  with  the  Fourier  series  computations  and  Major  Robinson  for  his  guidance  and 
encouragement  through  the  thesis  process.  Thanks  must  also  go  to  Dr.  Mark  Oxley 
for  his  assistance  with  the  mathematics.  I  am  especially  grateful  to  my  thesis  advisor, 

Dr.  Robert  Calico,  for  the  wisdom  and  insight  he  has  provided  me  through  the  past 
year.  Finally,  I  would  like  to  thank  my  wife,  Barbara  for  her  love  and  patience  during 
this  trying  time. 


James  Walter  Cole 
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Abstract 

The  stability  of  a  spinning  symmetric  satellite  in  an  elliptical  orbit  is  analyzed 
using  phase  planes  and  surface  of  section  techniques.  The  equations  of  motion  for 
the  satellite  attitude  are  presented  in  a  linearized  and  a  non-linear  form.  Floquet 
Theory  is  applied  to  the  development  of  a  control  system  for  two  unstable  modes  of 
the  satellite.  A  scalar  control  is  applied  using  angle  rates  as  feedback.  The  stability  of 
the  new  controlled  system  is  examined  in  controlled  modal  coordinates.  Comparisons 
of  the  linear  and  non-linear  system  motions  are  made  relative  to  changes  in  the 
control  gains.  Potential  chaotic  motion  limits  the  controller  gains  of  the  non-linear 
system.  ~T\  v,  <-  v,  ^  ,  '■ 
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Control  and  Stability  of  a  Spinning  Symmetric  Satellite  in  an 

Elliptical  Orbit 


I.  Introduction 

We  are  demanding  more  accurate  control  and  attitude  stability  from  our  satel¬ 
lites  while  sending  them  into  less  stable  environments.  The  control  design  developed 
here  is  relatively  easy  to  implement,  yet  it  compensates  for  a  complex,  dynamic  en¬ 
vironment.  This  fits  well  into  the  control  systems  for  satellites  going  into  elliptical 
orbits  or  in  other  periodic  systems  that  fit  the  form  of  the  Floquet  problem.  The 
non-linear  example  shows  that  some  caution  is  warranted  as  to  the  limits  of  con¬ 
trol  in  non-linear  systems.  In  some  configurations  non-linear  behavior  can  drive  the 
system  into  instability  or  chaos. 

The  stability  of  a  spinning  symmetric  satellite  was  examined  by  Kane  and 
Barba  [8:402-405].  The  development  of  the  equations  of  motion  from  their  work  is 
presented  in  Appendix  I.  They  applied  Floquet  Theory  to  determine  the  stability  of 
the  linearized  system.  Calico  and  Yeakel  developed  a  control  system  based  on  Flo¬ 
quet  Theory  to  control  one  unstable  mode  of  the  same  class  of  problem  [3:315-318]. 
Calico  and  Wiesel  extended  this  to  control  two  unstable  modes  [2:671-676],  Myers 
examined  several  types  of  controllers  for  this  kind  of  system  and  made  some  rough 
evaluation  of  the  performance  of  each  controller  type  [12:56-59].  Shell  examined  the 
scalar  controller  for  a  satellite  with  two  unstable  modes  using  phase  plane  techniques 
[  1 4:3.24-3.44] .  The  true  performance  limits  of  the  controller  are  still  unknown.  There 
are  no  guidelines  for  selecting  the  proper  parameter  values  for  this  controller  design. 
This  thesis  addresses  these  problems. 


The  theory  and  concepts  developed  here  can  be  applied  to  any  system  that 
fits  the  form  of  the  Floquet  problem.  The  example  problem  of  a  spinning  symmetric 
satellite  in  an  elliptic  orbit  is  well  defined  and  exhibits  some  interesting  characteris¬ 
tics.  A  full  description  of  the  coordinates  and  development  of  the  equations  of  motion 
are  given  in  Appendix  1.  The  equations  of  motion  can  be  expressed  in  matrix  form 
as 

x  —  A(t)x  (1.1) 

where  <4(r)  is  composed  of  periodic  elements  that  all  have  the  same  period.  The 
solution  to  the  this  system  of  equations  and  the  development  of  the  controller  involves 
the  application  of  Floquet  Theory. 
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II.  Theory 


This  chapter  examines  the  equations  of  motion  of  a  spinning  inertially  symmet¬ 
ric  satellite  in  an  elliptical  orbit  around  a  symmetrical  attracting  body,  determines 
the  satellite  stability,  develops  a  controller  and  examines  the  stability  of  the  con¬ 
trolled  system.  The  satellite  spin  axis  is  along  the  inertial  axis  of  symmetry  and 
perpendicular  to  the  orbit  plane,  Figure  2.1.  The  attitude  equations  of  the  satellite 
are  developed  and  linearized  about  an  equilibrium  solution  in  Appendix  1.  These 
linearized  attitude  equations  describe  a  linear  time  periodic  system.  Floquet  The¬ 
ory  is  applied  to  determine  the  stability  of  this  system.  A  scalar  controller  is  then 
developed  to  control  two  unstable  modes  of  the  linearized  system.  The  controller 
uses  two  rate  variables  from  the  system  state  vector  as  feedback  inputs.  Floquet 
Theory  is  then  applied  to  the  controlled  system  to  analyze  the  stability.  The  system 
is  evaluated  in  a  new  controlled  modal  coordinate  system,  derived  from  the  Floquet 
solution.  This  method  decouples  the  four  modes  of  the  system  and  provides  an  easy 
method  to  determine  system  stability. 

2. 1  Floquet  Theory 

The  solution  to  a  system  of  linear  time-periodic  equations  was  developed  by 
Floquet  in  the  late  1800’s.  Such  a  system  in  first  order  form  is  given  by 

x  =  ,4(r)x  (2.1 ) 

where  A(t)  =  A(t  +  T).  Any  set  of  independent  solutions, <p(r)  the  fundamental 
matrix,  must  satisfy  the  differential  equation 

<t>{ T )  =  A(r)d>(r)  (2.2) 
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Figure  2.1.  Orbit  Reference  Frame 
Floquet  found  that  the  fundamental  matrix  can  be  expressed  as 

<Kt)  =  P(r)err  (2.3) 

where  P  is  a  periodic  matrix  with  the  period  T  and  T  is  a  constant  matrix.  A  more 
convenient  and  equivalent  form  can  be  developed  through  a  similarity  transformation 
using  a  constant  non-singular  matrix  B  such  that  T  is  transformed  to  the  Jordan 
form, 

V»(r)  =  B~l<t>{T)B  =  F{r)ejT  (2.4) 

Then  rp  is  also  a  fundamental  matrix  and  F(r)  is  periodic,  with  a  period  of  T.  That 
is 

F(t)  =  F(t  +  T)  (2.5) 
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and  J  is  the  Jordan  form  of  T. 

Examine  equation  2.3  at  the  end  of  one  period, 


since  P(t)  is  periodic 


4>{t  +  T )  =  P(T  +  T)er<T+T> 


P{r )  =  P{t  +  T) 


(2.6) 


(2.7) 


substituting  this  in  to  the  equation  2.6 


4>(t  +  T)  =  P(r)erTerr  =  <p(r)e 


ST 


C,  a  constant  matrix,  replaces  erT  and  is  called  the  monodromy  matrix. 


4>{t  +  T)  =  4>(t  )err  =  <p(r)C 


The  monodromy  matrix,  C,  can  be  found  by  letting  r  =  0 


m~l<t>(T)  =  C 


If  o  is  the  principal  fundamental  matrix 


(2.8) 


(2.9) 


(2.10) 


4>(0)  =  I  (2.11) 

where  /  is  the  identity  matrix.  C  can  then  be  found  by  numerical  integration  of 
equation  2.2  using  the  initial  conditions  from  equation  2.11.  The  diagonal  elements 
of  J,  pi,  are  the  eigenvalues  of  T  and  are  called  characteristic  exponents.  The  system 
stability  can  be  determined  from  the  value  of  the  characteristic  exponents.  If  they 
have  negative  real  parts,  the  system  is  stable.  If  any  have  positive  real  parts,  the 
system  is  unstable.  The  J  matrix  can  be  determined  from  the  monodromy  matrix  C . 
The  eigenvalues  of  C,  Ai  (characteristic  multipliers),  are  related  to  the  eigenvalues 
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of  r  by 


Pi  =  J ln 


(2.12) 


This  analysis  will  allow  us  to  determine  the  system  stability  and  is  a  typical  applica¬ 
tion  of  Floquet  Theory.  But,  to  implement  a  controller  based  on  the  solution  to  the 
original  system  [2:671-676]  we  must  solve  for  the  F  matrix.  To  find  the  differential 
expression  for  F,  differentiate  equation  2.4  with  respect  to  r. 

0'(r)  =  F'{T)eJr  +  F{r)JejT  (2.13) 

We  note  that  every  fundamental  matrix  must  satisfy  the  differential  equation  2.2. 

*P'{  r )  =  A(t  )ip(r )  (2.14) 


Replacing  ip'  and  0  in  equation  2.14  with  the  expressions  from  equations  2.4  and 
2.13 

F'{r)eJr  +  F{r)JejT  =  A(r)F{T)eJr  (2.15) 

Dividing  both  sides  by  eJr  and  rearranging  results  in 


F'  =  AF-  FJ 


(2.16) 


To  find  the  initial  conditions  we  write 

i(r)  =  0(r)0(O)_1x(O)  (2.17) 

where  x(0)  is  the  initial  state  vector  and  ip  is  any  fundamental  matrix.  From  Floquet 
Theory  ,  any  fundamental  matrix  can  be  expressed  in  the  form  of 

0(r)  =  F(t )eJr  (2.18) 
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It  has  already  been  shown  how  to  find  the  J  matrix,  whose  diagonal  entries  are  the 
eigenvalues  of  the  C  matrix.  To  find  an  expression  for  the  F(0)  matrix,  take  the 
inverse  of  equation  2.18 

=  e~Jr  F(t)~1  (2.19) 

At  t  =  0 

4i0)~l  =  e°F(0)-1  =  F(  0)-'  (2.20) 

Then  use  this  value  in  equation  2.17 

x(r)  =  F(r)eJrF(0)-1x(0)  (2.21) 

Evaluated  at  t  =  T 

x(T)  =  F(T)eJTF(0)"1x(0)  (2.22) 

Writing  equation  2.17  in  terms  of  the  principal  fundamental  matrix 

x(r)  =  4>(t  )<^>(0)~1x(0)  (2.23) 

Evaluated  at  t  =  T 

x(T)  =  <t>{T)<f>{0)-lx{0)  (2.24) 

From  equation  2.10 

<t>{T)<t>{0)-1  =  <i>(T)I  =  4>(T)  =  C  (2.25) 

where  C  is  the  constant  monodromy  matrix.  Now  equate  equations  2.22  and  2.24 

F(T)e'/TF(0)-1x(0)  =  Cx(0)  (2.26) 

performing  the  necessary  linear  algebra  this  becomes 

F(T)eJT-CF{  0)  =  0  (2.27) 
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Noting  that  F(t)  is  periodic,  F[T)  =  F(0),  this  becomes 

F{ 0)eJT  -  CF(0)  =  0  (2.28) 

Since  J  is  diagonal,  this  can  be  written  as  a  set  of  i  equations  where  i  identifies  a 
column  in  the  F( 0)  matrix. 


9ifi( 0)  -  CM 0)  =  6  (2.29) 

Where  o,  is  the  ith  diagonal  element  of  eJT  matrix.  Or  expressed  as 

[gl  —  C]  /(0)  =  0  (2.30) 

which  is  an  eigenvalue  problem.  Each  column  of  the  F(0)  matrix  is  an  eigenvector 
of  the  monodromy  matrix  C.  This  will  provide  us  with  the  initial  conditions  for  the 
F  matrix.  As  previously  noted,  F  is  periodic,  and  if  numerically  integrated  over  one 
period,  a  set  of  Fourier  series  coefficients  can  be  extracted  from  the  data  set  and 
used  to  represent  the  elements  of  the  F  matrix.  Each  element  of  F  is  represented  as 
a  Fourier  series.  Additional  information  on  this  technique  can  be  found  in  reference 
[6:108-114]. 

Both  F  and  J  matrices  may  contain  complex  elements.  These  can  be  difficult 
to  work  with  but  equivalent  forms  can  be  developed  as  follows. 

•  F  matrix 

1.  If  the  column  vector,/,,  is  real  no  change  is  necessary. 

2.  If  a  pair  of  the  column  vectors  are  complex  conjugates,  separate  the  real 
and  imaginary  components  into  separate  vectors.  That  is 

Ireal+imagifreal—imag  ^  freah  fimag 

•  J  matrix 
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1.  If  the  eigenvalue,  p,,  is  real  no  change  is  necessary. 

2.  If  there  are  two  eigenvalues  that  form  a  complex  pair,  J  can  be  written 


with  the  diagonal  in  the  form  of  blocks 


Preal+imag 

0 

=> 

Prcal 

Pimag 

0 

Preal  —  imag 

Pimag 

preal 

These  new  equivalent  forms  of  F  and  J  contain  only  real  elements  and  are  much 
easier  to  deal  with  numerically. 

3.3  Control  Theory 

This  section  applies  Floquet  theory  to  develop  an  active  control  system  for  a 
spinning  symmetrical  satellite  in  an  elliptical  orbit.  The  equations  of  motion  for 
this  system  are  developed  in  Appendix  1.  Through  the  use  of  modal  coordinates, 
a  method  is  demonstrated  to  stabilize  any  system  that  has  a  two  unstable  modes. 
The  development  of  a  modal  control  system  for  this  system  differs  from  a  constant 
coefficient  linear  system  due  to  the  periodic  coefficients  of  the  differential  equations. 

The  first  step  in  developing  an  active  control  system  is  to  transform  the  physica1 
coordinates  to  modal  coordinates.  The  original  coordinates,  xT  =  |#i,  02,  0i,  ^2},  are 
related  to  the  modal  coordinates,  fj,  by  the  F  matrix  such  that 

x(r)  =  F(t)t}(t)  (2.31) 

Where  F(t)  is  the  periodic  matrix  that  solves  the  Floquet  problem.  Replacing  x  in 
equation  2.3  yields 

(F(t)t,(t))' =  A(t)F(t)t,(t)  (2.32) 

applying  the  chain  rule  to  the  derivative  and  reduction  yields 

n'(r)  =  F-'(t)  [A(t )F(t)  -  F\t)\  ,t(T)  (2.33) 
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From  equation  2.16,  replace  F' 


t)'{t )  =  F  1{t)  [ A(t)F(t )  -  A(t)F(t)  +  F(t)J}t]{t)  (2.34) 

This  reduces  to 

t/(t)  =  Jtj(t)  (2.35) 

The  eigenvector  matrix,  F(r),  transforms  the  time  periodic  system  to  a  constant- 
coefficient  system.  The  new  system  variables,  are  referred  to  as  modal  vari¬ 
ables  since  the  transformation  in  equation  2.31  performs  a  decoupling  similar  to  a 
constant-coefficent  system. 

The  controller  examined  here  was  developed  by  Calico  and  Wiesel  [2:673-674], 
It  is  a  two  mode  scalar  controller.  The  system  equation  with  a  generalized  controller 
ran  be  written  as 

x' =  A(t)x  +  Bu(t)  (2.36) 

Where  the  A  matrix  is  from  the  system  description  of  the  uncontrolled  system  and 
the  B  matrix  distributes  the  control  in  the  system.  This  can  also  be  expressed  in 
the  modal  variables. 

77'  =  Jfj  -f  F-1  Bu  (2.37) 

In  general,  B  has  the  fcrm 

0  •••  0 
0  •••  0 

l -.38) 

^32  •  •  •  63  n 

642  •  •  •  b4n 

The  number  of  columns  in  the  B  matrix  must  match  the  size  of  the  control  vector 
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u.  Each  element,  6(J,  may  be  time  periodic.  For  the  two  mode  scalar  controller 


(2.39) 


Which  variables  are  included  in  the  feedback  and  the  gain  associated  with  each 
variable  is  determined  bv  the  control  vector  u.  It  is  defined  as 


u  = 


(2.40) 


The  general  form  of  the  time  dependent  gain  matrix  is 


A'(r)  = 


A'u(r)  Kn(r)  K13(t)  A'm(t) 
A'2,(t)  Kn(r)  A'23(r)  K24(t) 

h'nl(T)  A'„2(t )  A'n3  ( T )  A'n4(r) 


(2.41) 


Again  for  the  two  mode  scalar  controller,  a  reduced  form  of  the  K(r)  matrix  is  used 


A'(r)  =  A',  A'2  0  0 


(2.42) 


where  both  non-zero  elements  are  periodic  functions  with  a  period  of  T .  It  is  assumed 
that  the  first  two  modes  are  unstable  and  need  to  be  controlled.  Expanding  the 
control  term,  for  our  two  mode  scalar  controller  results  in 


F~lBu  = 


III  f 1 2  f  13  fu 

/21  /22  fi3  f24 

/31  f32  fl3  /34 

/41  /42  f 43  f 44 


A',  I<2  0  0  f,  (2.43) 
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where  /,j  is  an  element  of  F  *.  Applying  matrix  algebra  reduces  this  equation  to 


F~lBu  = 


/l3  +  /l4 
/23  +  f24 
f33  +  /34 
/43  +  f 44 


and  let  A,  =  /, 3  +  /,4  it  becomes 


I<1  I<2  0  0  fj 


(2.44) 


F~lBn  = 


K1R1 

K2R1 

0 

0 

Ai  A2 

A  2  /?2 

0 

0 

A'./?3 

A  2  A3 

0 

0 

A 1 R4 

A'2A4 

0 

0 

Substitute  this  result  into  equation  2.37  and  expand  the  J  matrix 


ff  =  Jif+F-'BKif  = 


h\R\  +  P\  h2Ri  0  0 

A  j  /?2  A  2  7?2  +  /?2  0  0 
A 1  /?3  A  2  A3  P3  0 

A 1 A4  A  2  A4  0  p4 


(2.45) 


(2.46) 


This  is  also  a  Floquet  problem  and  the  techniques  from  Floquet  Theory  can  be 
applied  again  to  determine  the  stability  of  the  closed  loop  system  and  controller. 
This  Floquet  problem  can  also  be  expressed  in  physical  coordinates,  x.  Substituting 
into  equation  2.36  for  u 

x'  =  A(t)x  +  BK{T)fj  (2-47) 

and  replacing  f) 


This  reduces  to 


x'  =  A(t)x  BK(t)tj 

(2.47) 

x'  =  A{t)x  +  BK[t)F~xx 

(2.48) 

x'=  [a(t)  +  BK(t)F~'}x 

(2.49) 
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This  is  of  the  same  form  as  equation  2.2  where  the  new  periodic  matrix  includes 
the  control  elements  and  also  has  a  period  of  one  orbit.  The  solution  set  to  this 
Floquet  problem  will  be  of  the  same  form  and  can  be  expressed  as 

0(t)  =  F'(r)eJ'T  (2.50) 

where  F"  is  periodic  and  bounded,  J*  is  diagonal  and  constant.  The  stability  of  the 
controlled  linearized  system  is  determined  by  the  diagonal  elements  of  the  J*  matrix. 
If  the  real  part  is  negative  then  the  system  is  stable,  if  positive  then  the  system  is 
unstable.  A  new  set  of  modal  variables  for  the  controlled  system  can  be  defined  by 
the  transformation 

x  =  F'if  (2.51) 

The  controller  design  problem  is  to  find  the  proper  A' (r)  matrix  so  the  Poincar/’e 
exponents  are  all  negative  [2:674]  [12:27],  Working  with  modal  variables  this  prob¬ 
lem  reduces  to  two  coupled  differential  equations. 

rh  =  ( I<iRi  +  <rK  +  (A'afli  +  "W  ('2-52) 

rh  ~  (XiR2  —  +  ( A 2 R2  +  v) V2  (2.53) 

where  px  and  p2  are  a  complex  pair  that  have  been  transformed  into  two  real  parts 
a  and  u>.  There  is  no  known  conventional  solution  method  for  these  two  equations. 
Since  the  system  is  periodic,  it  can  be  written  in  the  form  of  a  Floquet  problem. 

V  =  X(r)i j  (2.54) 


And  applying  another  theorem 


(2.55) 
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shows  a  relationship  that  provides  some  insight  into  the  correct  values  for  the  A' 
matrix.  Here  A,  are  the  characteristic  multipliers  and  /r[.Y(s)]  is  the  trace  of  the  .Y 
matrix.  If  the  constant  terms  in  the  Fourier  series  form  of  the  F  matrix  are  large 
enough  to  be  of  practical  use  in  the  feedback,  then  K\  and  A'2  can  be  constants,  and 
the  trace  becomes 


tr  [.Y(s)]  =  (A>i  +  K2a2  +  2a)  +  P(t)  (2.56) 

where  the  a,  are  the  sum  of  the  respective  constant  Fourier  terms  and  P(r)  is  the 
remaining  periodic  terms.  Now 

ln(AiA2)  =  jf  [Kicti  +  K2a2  +  2cr  +  P(t)]  dr  (2.57) 

Completing  the  integration  over  one  period  yields 

ln(AjA2)  =  ( Aiai  +  A2o;2  4-  2  a)T  (2.58) 

or 

—  ln(Aj )  +  —  ln(A2)  =  h\ot\  +  I\2a2  +  2a  (2.59) 

Note  that  ( 1/T*)  In  A,-  =  the  characteristic  exponent.  This  results  in 

P\  d"  P2  =  AjOfi  +  A2a2  +  2a  (2.60) 

For  the  characteristic  exponents  to  be  negative,  the  sum  must  also  be  negative. 
Unfortunately  the  converse  is  not  true.  We  cannot  guarantee  both  exponents  will  be 
negative  if  the  sum  is  negative,  except  in  the  case  of  a  complex  pair.  This  relationship 
will  define  a  starting  region  for  values  of  I\  but  cannot  always  give  stable  solutions. 

In  summary,  we  have  examined  the  stability  of  the  uncontrolled  system  using 
Floquet  Theory,  developed  a  controller  that  can  control  two  unstable  modes  with 
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the  proper  selection  of  the  feedback  gains,  and  also  looked  at  the  controlled  system 
stability  again  with  Floquet  Theory.  This  stability  determination  applies  only  to 
the  linearized  system  and  the  controller  can  only  be  shown  mathematically  to  be 
effective  controlling  the  linearized  system.  What  happens  when  it  is  applied  to  the 
non-linear  system?  This  will  require  some  special  analysis  techniques. 
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III.  Stability  Analysis 


This  chapter  describes  the  analysis  techniques  used  to  examine  the  linearized 
and  the  non-linear  system.  Some  of  the  tools  have  been  known  since  the  late  1800’s 
but  one  of  the  most  powerful  tools  is  a  recent  development,  the  digital  computer. 
For  without  digital  computer  and  graphical  output  devices  the  other  tools  are  almost 
unusable.  Even  in  relatively  simple  problems,  like  our  example  with  four  differential 
equations,  computers  and  graphical  output  is  essential  to  compute,  manipulate  and 
display  the  data. 

3. 1  Classical  Methods 

In  any  undergraduate  text  on  control  methods  for  linear  systems  a  few  basic 
methods  of  determining  the  stability  of  a  system  will  be  presented.  These  are  quite 
effective  for  linear  system  design  and  performance  evaluation.  Two  principal  methods 
that  apply  graphic  output  are  the  root  locus  and  Bode  plots. 

3.1.1  Root  Locus  The  root  locus  plot  describes  the  path  of  the  roots  in  the 
complex  plane  as  some  parameter,  usually  the  feedback  gain,  is  varied.  This  is  an 
effective  design  tool  for  linear  systems.  Identifying  the  source  of  instability  and  the 
selection  of  a  gain  that  will  stabilize  the  system  is  fairly  easy  for  linear  systems.  With 
the  use  of  Floquet  theory  and  modal  coordinates,  this  technique  can  be  applied  to 
the  example  problem  by  plotting  the  characteristic  exponents  on  the  complex  plane 
as  the  gain  of  the  feedback  circuit  is  changed. 

3.1.2  Bode  Plots  Bode  plots  provide  a  display  of  the  frequency  response  of 
a  linear  system.  The  gain  and  phase  of  the  system  are  plotted  versus  the  input 
frequency.  Complete  design  and  stability  can  be  determined  from  the  two  plots. 
But  as  with  the  root  locus,  it  too  is  limited  to  linear  constant  coefficient  systems. 
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Modal  coordinates  can  be  applied  to  linear  time  periodic  systems.  Through  a 
modal  coordinate  transformation  the  periodic  system  can  be  represented  as  a  linear 
constant  coefficient  system.  This  allows  all  of  the  linear  methods  to  be  applied  to  a 
linear  time  periodic  system.  Floquet  theory  provides  that  type  of  transformation. 

3.2  Son-linear  Methods 

Several  methods  have  been  developed  that  can  be  applied  to  non-linear  sys¬ 
tems.  Often  they  do  not  provide  detailed  design  information  but  they  show  trends 
in  the  behavior  and  define  the  regions  of  stability.  These  techniques  can  be  applied 
to  linear  systems  but  the  specialized  linear  methods  give  precise  design  information. 

3.2.1  Phase  Space  The  use  of  phase  space  to  examine  non-linear  systems 
provides  a  geometric  representation  of  the  system  moti^-  ’"lach  point  of  the  phase 
space  represents  a  system  state.  A  time  trace  in  the  phase  space  shows  the  trans¬ 
formation  of  the  system  from  one  state  to  another.  A  typical  example  of  a  two 
dimensional  phase  space  is  a  pendulum,  where  the  position  (0),  and  velocity  (0),  are 
used  as  a  Cartesian  coordinate  system,  Figure  3.1.  This  provides  a  two  dimensional 
space  or  phase  plane  to  examine  the  motion.  The  origin  of  the  axis  represents  the 
equilibrium  position,  the  pendulum  is  stationary.  It  is  surrounded  with  a  series  of 
concentric  rings,  each  representing  a  stable  periodic  motion.  As  one  examines  the 
phase  plane  further  from  the  equilibrium  point,  the  non-linearity  of  the  system  be¬ 
comes  apparent.  Other  equilibrium  points  are  located  at  values  of  0  that  correspond 
to  a  vertical  position  of  the  pendulum,  these  are  all  unstable.  If  the  pendulum  is 
given  enough  energy  it  will  swing  completely  around  the  pivot  point.  This  motion 
is  represented  by  the  wavy  lines  outside  of  the  closed  orbits. 

This  is  a  simple  non-linear  system  with  no  damping  or  forcing  input.  It  can 
easily  be  shown  in  only  two  dimensions,  and  yet  the  phase  plane  shows  complex 
motion.  Higher  dimensional  systems  with  damping  and  forcing  functions  can  exhibit 
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highly  complex  motion  in  their  respective  phase  space. 

On  the  phase  plane  there  several  geometric  patterns  that  can  be  classified 
according  to  the  behavior  of  the  trajectories  in  the  vicinity,  Figure  3.2.  These 
structures  represent  solutions  to  a  two  dimensional  system. 

•  A  node  is  characterized  on  the  phase  plane  by  the  intersection  of  trajectories 
from  the  surrounding  space.  They  may  be  either  stable  or  unstable.  If  the 
trajectories  move  toward  the  node  as  time  progresses  then  it  is  a  stable  node. 
If  the  trajectories  move  away  from  the  node  as  time  moves  forward,  the  node 
is  considered  unstable. 

•  A  center  is  an  equilibrium  point  surrounded  by  closed  path  trajectories. 

•  An  orbit  surrounds  a  center  and  represents  a  periodic  solution. 

•  A  saddlepoint  is  an  unstable  solution  that  is  identified  by  the  intersection  of  two 
or  more  trajectories.  The  distinguishing  characteristic  is  half  of  the  trajectories 
approach  the  solution  and  the  other  half  of  the  trajectories  move  away  from 
the  point. 
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•  A  spiral  is  a  special  type  on  node.  The  trajectories  do  not  actually  intersect 
but  spiral  around  the  node.  Spirals  can  be  stable  or  unstable. 

•  A  limit  cycle  is  a  closed  trajectory  and  represents  a  periodic  solution.  It  differs 
from  an  orbit  in  that  the  surrounding  trajectories  will  converge  or  diverge  to 
the  limit  cycle  orbit.  If  they  converge  to  the  limit  cycle  it  is  considered  a  stable 
limit  cycle.  If  the  trajectories  diverge  from  the  orbit,  it  is  unstable. 

Stable  structures  are  also  known  as  attractors  and  the  unstable  structures  are 
repellers  because  of  the  way  the  phase  paths  are  attracted  or  repelled  from  the  regions 
of  the  structure.  All  of  these  structures  can  appear  on  the  phase  plane  in  various 
combinations.  For  a  full  discussion  of  these  structures  see  references  [10:170-206] 
and  [7:1-54], 

A  domain  of  attraction  can  be  defined  for  each  stable  structure  [15:49].  This 
is  the  region  of  the  phase  space  that  trajectories  lead  to  that  stable  solution.  If  all 
trajectories  in  the  phase  space  lead  to  a  solution  then  the  domain  of  attraction  is 
the  entire  phase  space  and  it  is  considered  globally  stable. 

Stable  structures  are  easily  identified  by  the  converging  flow  of  the  trajectories 
in  the  surrounding  phase  space.  Unstable  structures  can  be  identified  by  examining 
the  flow  of  trajectories  in  a  negative  time  period.  This  is  usually  accomplished  by 
running  the  numerical  integrator  with  a  negative  time  step.  The  unstable  nodes, 
spirals  and  limit  cycles  can  be  easily  identified  and  appear  as  stable  structures  in 
reverse  time.  The  behavior  of  an  orbit  or  a  center  will  be  unchanged  and  a  saddle 
point  will  still  be  unstable  and  exhibit  a  complementary  structure. 

Phase  space  is  a  useful  tool,  especially  when  the  system  can  be  represented  with 
only  two  variables.  Higher  dimensional  phase  space  is  difficult  to  model  and  visualize, 
however  there  are  other  methods  that  will  provide  insight  into  higher  dimensional 
phase  space. 
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Stable  Limit  Cycle  Unstable  Limit  Cycle 


Figure  3.2.  Example  Phase  Plane  Structures 
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3.2.2  Poincare  Maps  A  Poincare  map,  also  referred  to  as  a  Poincare  section, 
is  constructed  by  defining  a  plane  in  three  or  higher  dimensional  phase  space.  Every 
time  the  trajectory  passes  through  the  plane,  it  is  mapped  as  a  point  on  the  plane. 
If  many  trajectories  are  mapped  in  this  way,  the  underlying  behavior  of  the  system 
can  be  determined.  According  to  Dr.  Francis  C.  Moon,  the  Poincare  maps  can  be 
categorized  by  the  patterns  that  are  generated  [11:53].  These  structures  are  observed 
on  a  two  dimensional  plane  in  the  phase  space  that  the  trajectories  are  allowed  to 
puncture  as  they  travel  in  a  three  or  higher  dimensional  space. 

•  A  finite  number  of  points  indicates  a  periodic  motion  or  harmonic  oscillation. 

•  Closed  curves  indicate  a  quasiperiodic  motion  or  a  combination  of  two  periodic 
motions  that  are  not  harmonic. 

•  A  fractal  shaped  curve  is  a  strange  attractor  in  a  higher  dimension  phase  space. 

•  A  fuzzy  or  dispersed  collection  of  points  has  several  possible  causes  due  to 
the  lack  of  definition.  It  could  be  a  system  with  three  or  more  non-harmonic 
frequencies,  a  lightly  damped  strange  attractor  or  a  strange  attractor  in  more 
that  three  dimensions.  Other  tests  will  have  to  be  used  to  determine  the  true 
nature  of  the  system. 

3.2.2. 1  Damping  Effects  A  Hamiltonian  system  is  conservative  and  has 
no  damping.  The  Poincare  maps  for  Hamiltonian  systems  will  tend  to  fill  a  region  on 
the  phase  plane.  The  structure  of  the  attractor  will  be  so  disperse  that  no  pattern 
will  be  evident.  The  trend  is  for  damping  to  reduce  the  dispersion  of  the  map  and  if 
increased  sufficiently  it  will  cause  a  fractal  shape  to  reduce  to  a  single  curve  [11:63]. 

3. 2. 2. 2  Torus  Mapping  A  system  that  has  a  dominant  frequency  can 
be  mapped  on  the  surface  of  a  torus.  One  cycle  of  the  forcing  function  is  mapped 
as  a  polar  coordinate  from  the  center  of  the  torus.  The  phase  plane  of  the  Poincare 
map  is  perpendicular  to  the  torus  and  cuts  a  cross-section,  Figure  3.3.  The  motion 
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of  a  periodic  system  will  make  lines  around  the  torus,  similar  to  a  barber  pole.  A 
quasiperiodic  system  will  completely  cover  the  surface  of  the  torus  with  paths.  A 
chaotic  system  will  have  a  highly  deformed  torus  with  a  fractal  cross-section  on  the 
Poincare  section.  This  type  of  map  can  be  used  to  examine  strange  attractors  and 
quasiperiodic  attractors  as  they  change  in  relation  to  the  phase  of  the  dominant 
period.  By  taking  Poincare  maps  at  several  phases  in  the  cycle  of  the  dominant 
frequency  a  time  picture  of  the  strange  attractor  can  be  examined.  Just  move  the 
cross-section  of  the  Poincare  map  around  the  torus  to  create  new  maps. 

All  of  the  previous  techniques  apply  to  third  order  systems.  For  example,  the 
phase  space  coordinates  would  be  x,  x,  and  mod(u>t,  2n).  Dr.  Moon  has  demon¬ 
strated  a  method  referred  to  as  a  ’Double  Poincare  Map’  [IT.  142].  This  method  will 
find  strange  attractors  in  higher  dimension  phase  space  by  taking  multiple  cuts  in 
to  higher  dimension  torus.  He  demonstrates  the  method  on  a  fourth  order  system 
and  states  it  can  be  generalized  to  higher  order  systems. 

3.2.3  Power  Spectrum  The  power  spectrum  or  Fourier  spectrum  can  be  ap¬ 
plied  to  differentiate  between  a  periodic  system  and  a  chaotic  one.  The  power  spec¬ 
trum  will  identify  any  periodic  component  of  the  system  as  a  spike.  A  quasiperiodic 
system  will  have  multiple  spikes  at  the  basic  frequencies  involved  and  mixed  har¬ 
monics  of  the  basic  frequencies.  Chaotic  systems  will  have  a  broad  spectrum,  the 
regions  between  spikes,  if  any,  will  be  filled  [1:314-337], 

3.3  Chaos 

To  quote  James  Gleick,  ’Where  Chaos  begins,  classical  science  stops.’  [5:3]. 
Until  recently,  random  motion  was  an  area  that  had  to  be  dealt  with  stochastically.  A 
continuum  of  complexity  is  presented  by  Heinz  Pagels,  leading  from  highly  ordered  to 
seemingly  random  systems.  The  crystal  lattice  of  a  diamond  is  highly  structured  and 
predictable,  a  rose  presents  a  complex  structure  with  order  but  not  as  predictable. 


3-7 


3-8 


the  motion  of  gas  molecules  in  a  volume  presents  a  highly  complex  system  will  very 
little  order  and  difficult  to  predict  the  future  motion. 

In  these  highly  random  systems  one  can  speak  of  the  average  behavior  by 
the  use  of  statistics  [13:69-70],  Why  should  one  have  to  resort  to  statistics  in 
a  deterministic  problem?  For  a  dynamic  system  whose  motion  is  described  by  a 
system  of  ordinary  differential  equations,  given  a  set  of  initial  conditions,  the  future 
motion  of  the  system  is  prescribed.  There  is  nothing  left  of  chance. 

For  example,  one  of  the  simplest  of  dynamic  systems  that  is  considered  a 
"random  event”  is  the  result  of  a  coin  toss.  The  result  of  the  coin  toss  is  determined 
by  it’s  initial  velocity  and  spin  rate.  The  apparent  randomness  comes  from  ’sensitive 
dependence  on  initial  conditions’.  From  various  readings  this  is  the  basic  requirement 
for  a  chaotic  system.  Any  slight  variation  in  either  velocity  or  spin  rate  can  cause  a 
change  in  the  result  of  the  coin  toss. 

To  create  a  plot  of  phase  space  for  the  coin,  one  axis  is  the  velocity  and  the 
other  the  spin  rate.  Color  the  regions  that  result  in  heads  black  and  in  tails  white. 
There  will  be  areas  of  white  and  areas  of  black,  but  the  majority  of  the  phase  space 
will  consist  of  finely  intermixed  regions.  It  is  this  mixed  region  that  causes  us  to 
think  of  the  deterministic  system  as  being  random.  The  initial  conditions  cannot 
be  set  precisely  enough  to  be  able  to  predict  the  result.  The  small  changes  in  the 
initial  conditions  end  at  a  different  solutions  [13:64-68].  A  extreme  example  of  this 
is  found  in  trying  to  predict  the  weather.  Referred  to  as  ’The  Butterfly  Effect',  the 
system  is  so  sensitive  to  disturbance  a  butterfly  perturbing  the  air  in  Peking  today 
causes  storms  in  New  York  next  month  [5:20-23]. 

3.3.1  Strange  Attractors  Apparently  random  motion  from  a  deterministic 
system  was  studied  by  Lorenz  from  a  set  of  three  differential  equations  he  was  using 
as  a  weather  model.  His  discovery  of  the  strange  attractor  in  phase  space  lead  to 
new  understanding  of  deterministic  chaotic  systems.  This  started  a  revolution  in  the 
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analysis  of  non-linear  systems. 

An  attractor  is  a  region  in  phase  space  that  is  occupied  by  a  steady  state 
trajectory.  The  simplest  attractor  is  a  node,  with  dimension  of  0.  The  limit  cycle  is 
an  attractor  with  a  dimension  of  1.  A  quasiperiodic  torus  is  also  an  attractor  and 
has  a  dimension  of  2.  A  strange  attractor  is  a  region  or  surface  in  the  phase  space 
into  which  the  steady  state  solution  resides  and  has  a  non-integer  dimension.  The 
basin  of  attraction  fills  the  region  around  an  attractor.  Any  starting  point  in  the 
basin  of  attraction  will  result  in  a  steady  state  solution  to  that  attractor. 

A  strange  attractor  will  confine  the  solution  to  a  bounded  region  of  the  phase 
space.  The  trajectory  will  move  through  this  region  on  a  continuous  path  that  never 
intersects.  Any  two  trajectories,  arbitrarily  close,  in  this  region  will  diverge  rapidly 
over  time.  The  other  attractors  do  not  exhibit  this  property,  in  their  basins  of 
attraction  two  trajectories,  arbitrarily  close,  will  remain  close.  When  viewed  on  a 
Poincare  section  a  strange  attractor  will  have  a  fractal  shape  [4:200]. 

Transitions  to  Chaos  As  mentioned  earlier,  there  is  a  continuum  that 
exits  from  the  simple  to  the  complex.  In  this  realm  of  non-linear  dynamics,  there 
are  paths  that  a  system  can  take  in  transitioning  from  a  simple  stable  system  into  a 
chaotic  one.  Many  dynamic  systems  exhibit  period  doubling  or  bifurcations  as  the 
control  parameter  is  varied.  This  is  required  in  one-dimensional  systems  but  cannot 
be  generalized  as  a  requirement  for  higher  order  systems  but  is  sometimes  observed 
[15:126-128).  Higher  order  systems  have  been  observed  to  have  homoclinic  orbits. 
Homoclinic  orbits  have  the  property  that  from  an  initial  circular  starting  region,  as 
the  system  moves  in  time  the  initial  region  is  mapped  into  a  horseshoe  shape.  As 
time  continues  the  mapping  will  be  repeated.  This  results  in  a  finely  dispersed  map 
with  a  fractal  pattern  [11:161]. 
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IV.  Results 


In  this  chapter,  some  of  the  previous  work  on  the  example  problem  is  examined 
and  several  of  the  techniques  discussed  in  chapter  3  are  applied  to  reveal  new  stability 
information.  The  software  used  to  run  the  simulations  was  written  in  Fortran77  and 
Microsoft  Fortran  4.1  used  the  IMSL  math  libraries.  All  calculations  were  performed 
using  double  precision  routines.  Some  of  the  code  was  developed  by  Captain  Gregory 
Myers  and  then  expanded  by  Captain  Dale  Shell  to  include  the  non-linear  equations 
of  motion.  I  have  added  the  code  to  handle  real  as  well  as  complex  roots,  deal 
with  the  period  doubling  of  the  system  and  perform  the  controlled  modal  coordinate 
transformation.  The  simulations  of  the  linearized  system  were  run  on  an  Zenith  248 
microcomputer  under  MSDOS  operating  system.  The  non-linear  runs  were  executed 
on  an  Elxsi  6400  superminicomputer  under  the  EMBOS  operating  system. 

4-1  Uncontrolled  System 

The  uncontrolled  system  is  classified  as  a  Hamiltonian  system.  There  is  a 
conservation  of  energy  in  the  system.  It’s  behavior  was  examined  by  both  Myers 
and  Shell  in  both  physical  and  modal  coordinates. 

4- LI  Physical  Coordinates  Myers  found  that  the  linearized  equations  of  mo¬ 
tion  when  examined  in  physical  coordinates  exhibited  some  instability.  This  was 
determined  by  examining  the  parameter  <f>,  the  angle  between  the  satellite  spin  axis 
and  the  Az  coordinate  axis.  Over  a  period  of  six  orbits  the  size  of  <j>  grew  from  the 
perturbed  value  to  a  larger  average  value  as  shown  on  figure  4.1. 
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4-1-2  Modal  Coordinates  When  Myers  plotted  this  data  in  the  modal  coor¬ 
dinates  defined  by  the  Floquet  solution  as 

fj  =  F~x(t)  x  (4.1) 

Over  the  first  six  periods,  modes  one  and  two  appear  to  be  unstable  and  modes  three 
and  four  are  stable.  See  figures  4.2,  4.3,  4.4  and  4.5.  Shell  extended  this  work 
by  examining  the  non-linear  equations  of  motion  in  the  same  coordinate  system  and 
extending  the  integration  of  the  equations  to  twelve  orbits.  Still  it  appears  that 
the  uncontrolled  system  has  two  unstable  modes  and  two  stable.  He  did  note  the 
irregular  motion  of  the  stable  modes  as  they  approached  the  center  and  did  correctly 
identify  it  as  coupling  of  the  non-linear  terms  with  the  other  two  modes.  See  figures 
4.6  and  4.7. 

If  he  had  continued  the  integration  a  little  further  in  time,  a  undiscovered 
feature  of  the  system  would  have  presented  itself.  Figures  4.8  and  4.9  show  the 
uncontrolled  system  for  40  orbits.  This  is  more  than  sufficient  time  to  observe  full 
motion  of  the  system.  With  the  modes  lightly  coupled  there  is  an  energy  exchange 
that  occurs  between  the  modes.  Remember  this  is  a  conservative  system.  If  one  set 
of  modes  is  declared  stable  and  the  other  two  unstable,  the  energy  from  the  stable 
set  is  being  transferred  to  the  unstable  set.  As  the  t)\  and  772  modes  spiral  out,  the 
773  and  774  modes  spiral  inward  to  the  center.  The  non-linear  terms  are  apparent  in 
the  773  and  774  modes  first.  When  the  non-linear  terms  become  dominant,  it  causes  a 
reversal  in  the  energy  exchange.  r)X  and  772  make  a  very  irregular  path  back  to  a  low 
energy  state  while  773  and  774  move  along  an  irregular  path  to  a  high  energy  state. 
This  energy  exchange  continues  through  all  time  and  never  runs  down  or  becomes 
more  unstable  as  long  as  there  are  no  other  external  influences. 

4-1.3  Poincare  Map  As  stated  in  Chapter  3,  a  Poincare  map  of  a  conservative 
system  will  cover  a  region  of  phase  space  with  a  disperse  scatter  of  points.  See 
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figures  4.10  and  4.11.  This  a  new  technique  in  the  use  of  Poincare  maps,  the 
four  dimensional  system  is  broken  into  two,  two  dimensional  maps.  This  is  possible 
due  to  the  modes  being  independent  and  spanning  the  space  completely.  These  two 
maps  give  little  further  information  than  confirming  what  was  already  know  of  the 
non-linear  system. 

4-2  Controlled  System 

A  control  system  added  to  the  satellite  is  supposed  to  damp  out  the  pertur¬ 
bations  and  return  the  system  to  the  preferred  equilibrium  position.  Examining  the 
controlled  system,  we  have  two  goals 

•  determine  how  effective  the  control  system  is  or  how  large  of  a  perturbation 
will  it  correct 

•  determine  the  proper  setting  of  the  controller  gain  and  examine  how  it  influ¬ 
ences  the  performance  of  the  controller 

4-2.1  Physical  Coordinates  Myers  examined  the  linearized  system  in  his  the¬ 
sis  and  was  able  to  determine  a  method  for  selecting  a  gain  level  for  a  stable  system 
and  projected  that  a  ’’most  stable"  solution  would  be  found  at  the  point  where  the 
complex  roots  intersect  with  the  real  axis.  He  also  successfully  selected  gain  levels 
that  stabilized  the  linearized  system  and  plotted  the  system  response  in  physical 
coordinates  [12:56].  The  use  of  physical  coordinates  shows  little  detail  of  the  system 
behavior  other  that  an  general  trend  of  the  decrease  in  magnitude  of  the  perturbed 
motion. 

4-2.2  Root  Locus  Myers  also  sketched  a  root  locus  plot  of  the  control  systems 
he  examined  but  he  never  actually  computed  the  necessary  data  to  make  an  exact 
plot.  There  is  a  hidden  difficulty  in  this  system,  the  phenomena  know  as  period 
doubling  occurs  in  this  periodic  system.  This  can  be  resolved  without  too  much 
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difficulty  if  it  is  noted  that  any  system  that  has  a  period  of  2k  can  be  treated  as  if 
it  had  a  period  of  47T.  Figure  4.12  is  a  plot  of  the  Poincare  exponents  as  a  function 
of  gain  for  the  linearized  system.  The  transition  from  a  27r  period  to  4k  occurs  at 
the  intersection  of  the  real  axis.  I  should  point  out  that  there  is  one  false  point  on 
the  plot  at  the  intersection  of  the  real  axis.  I  had  to  ’create’  this  one  to  make  the 
plotting  package  cause  an  intersection  of  the  two  complex  roots  with  the  axis. 

A  close  examination  shows  that  with  a  gain  setting  of  zero,  there  are  two 
stable  roots  and  two  unstable  roots,  as  the  gain  is  increased  on  the  controller  the 
two  unstable  roots  cross  the  imaginary  axis.  When  the  gain  is  about  0.2  the  roots 
intersect  at  the  real  axis,  at  0.4  there  are  two  distinct  real  roots. 

This  does  provide  some  detailed  information  about  the  linearized  system.  Clas¬ 
sical  control  techniques  can  be  used  to  optimize  the  system  for  the  desired  perfor¬ 
mance  based  on  this  kind  of  plot.  This  will  help  us  select  the  gain  required  for  our 
controller  but  it  falls  short  in  predicting  controller  performance  in  the  non-linear 
system.  As  long  as  the  system  is  restricted  to  the  region  around  the  equilibrium 
point,  it  can  be  expected  to  perform  in  a  near  linear  manner. 

4.2.3  Modal  Coordinates  Shell,  in  his  examination,  used  the  modal  coordi¬ 
nates  and  phase  planes  to  examine  the  system  response.  He  examined  the  linearized 
and  the  non-linear  equations  of  motion  in  the  modal  phase  plane.  He  limited  his 
examination  to  a  fixed  gain  level  and  just  varied  the  perturbation  to  the  system  to 
determine  its  response.  His  use  of  the  modal  coordinates  for  the  uncontrolled  system, 
defined  by 

f)  =  F~x{t)x  (4.2) 

was  not  effective  in  decoupling  the  modes  of  the  controlled  system.  This  requires 
the  use  of  a  controlled  modal  coordinate  system  defined  as 

rf  =  F*~1(t)x  (4.3) 
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Figure  4.12  Root  Locui  of  Linearized 


Table  4.1.  Starting  values  for  linearized  system 


m 

V3 

T)4 

A 

0.5 

0.0 

0.0 

0.0 

B 

1.0 

0.0 

0.0 

0.0 

C 

1.5 

0.0 

0.0 

0.0 

D 

2.0 

0.0 

0.0 

0.0 

E 

3.0 

0.0 

0.0 

0.0 

F 

4.0 

0.0 

0.0 

0.0 

G 

5.0 

0.0 

0.0 

0.0 

H 

6.0 

0.0 

0.0 

0.0 

I 

7.0 

0.0 

0.0 

0.0 

that  can  be  derived  from  the  solution  to  the  Floquet  problem  presented  by  the  con¬ 
trolled  system,  equation  2.49.  This  will  define  a  system  of  independent  coordinates 
that  span  the  phase  space  and  decouple  the  modes  of  the  linearized  controlled  system. 
It  should  also  be  noted  that  it  is  invertible  and  there  is  a  one  to  one  correspondence 
between  the  physical  coordinates,  x,  and  the  modal  coordinates,  fjm. 

If  the  gain  is  set  at  0.4,  as  Shell  did  in  his  analysis,  the  system  performance 
can  be  examined  in  the  phase  space  of  the  controlled  modal  coordinates.  First  the 
linearized  system  is  examined  using  the  starting  values  in  Table  4.1.  The  phase  plane 
traces  are  shown  in  Figures  4.13  and  4.14.  These  figures  present  several  interesting 
features  of  the  linearized  system.  The  controller  works  over  the  entire  phase  space  in 
a  predictable  manner.  All  four  modes  are  stable  and  in  a  linear  system  the  designer 
could  proceed  to  optimize  the  gain  for  the  desired  performance  features.  Also  note 
that  the  equilibrium  point  is  transformed  from  x  =  [0,0,0,0]T  to  fj  =  [0,0,0,0]T. 
This  is  not  true  for  other  points  in  the  two  spaces. 

This  gives  very  little  additional  information  than  the  root  locus.  It  was  already 
known  the  linearized  system  would  be  stable  for  this  value  of  the  gain  and  the  phase 
portraits  of  the  system  support  that  fact.  Let  us  now  turn  to  the  non-linear  equations 
of  motion  and  examine  them  in  the  phase  space.  To  start,  select  a  perturbation  point 
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Table  4.2.  Starting  values  for  0  quadrants 


Ox 

02 

n 

mm 

A 

0.15 

0.15 

0.0 

0.0 

B 

-0.15 

0.15 

0.0 

0.0 

C 

0.15 

-0.15 

0.0 

0.0 

D 

-0.15 

-0.15 

0.0 

0.0 

Table  4.3.  Starting  values  for  r]\  axis 


0i 

02 

03 

04 

A 

0.5 

0.0 

0.0 

0.0 

B 

0.6 

0.0 

0.0 

0.0 

C 

0.7 

0.0 

0.0 

0.0 

D 

0.8 

0.0 

0.0 

0.0 

E 

0.9 

0.0 

0.0 

0.0 

F 

.0 

0.0 

0.0 

0.0 

G 

1.1 

0.0 

0.0 

0.0 

H 

1.2 

0.0 

0.0 

0.0 

in  each  quadrant  of  Oi  and  02  in  the  physical  phase  space.  The  starting  values  are 
presented  in  Table  4.2  and  the  phase  portraits  are  on  Figures  4.15  and  1.16. 
These  paint  a  different  picture  from  the  linearized  system.  As  expected,  in  the 
region  of  the  equilibrium  point  the  system  shows  behavior  almost  identical  to  the 
linearized  system,  but  in  the  regions  further  out  it  is  clearly  non-linear  and  does  not 
immediately  show  any  signs  of  being  stable  and  may  even  be  chaotic. 

Exploring  the  phase  space  at  several  points  along  the  r;,  axis,  at  about  i]  = 
[l,0,0,0]r  the  traces  no  longer  appear  to  converge  to  the  equilibrium  point,  Table 
4.3  and  Figures  4.17  and  4.18.  A  similar  adventure  along  the  t]2  axis  shows  a  similar 
state  is  reached  at  about  17  =  [0, 1,0, 0]r,  Table  4.4  and  Figures  4.19  and  4.20. 

The  operational  limits  of  the  controller  are  slowly  being  outlined  but  this  is  not 
the  final  answer,  a  period  of  12  orbits  is  not  sufficient  to  determine  the  final  trajectory 


4-21 


Figure  4.15  Non-iinesr  ij,  ij,  for  12  orbits 
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Figure  4.19  Non-lineer  jj.  ve  v2  tor  12  o^bite 


Table  4.4.  Starting  values  for  t?2  axis 


Vi 

n-i 

V3 

V4 

A 

0.0 

0.5 

0.0 

0.0 

B 

0.0 

0.6 

0.0 

0.0 

C 

0.0 

0.7 

0.0 

0.0 

D 

0.0 

0.8 

0.0 

0.0 

E 

0.0 

0.9 

0.0 

0.0 

F 

0.0 

1.0 

0.0 

0.0 

G 

0.0 

1.1 

0.0 

0.0 

H 

0.0 

1.3 

0.0 

0.0 

of  the  points  started  in  the  non-linear  region.  Start  a  trajectory  at  fj  =  [0.9, 0, 0, 0]r 
and  integrate  it  for  200  orbits,  the  results  are  shown  in  Figures  4.21  and  4.22. 
These  figures  indicate  the  controller  is  effective  at  this  condition. 

At  a  starting  point  of  fj  =  [0.95,0,0,0]r  the  results  are  completely  different  as 
shown  in  Figures  4.23  and  4.24.  The  first  few  orbits  are  found  on  the  interior  of 
this  structure  and  then  the  system  appears  to  settle  into  an  orbital  type  structure. 
This  would  r>~t  be  desirable  for  a  satellite  control  system,  it  would  quickly  exhaust 
its  fuel  supply. 

There  is  a  region  around  the  equilibrium  point  that  exhibits  a  stable  behavior 
and  if  it  can  be  assured  of  not  exceeding  those  bounds,  the  controller  will  perform  its 
intended  function.  Examining  points  outside  of  this  orbital  structure,  for  example 
f)  =  [3,0,0, 0]r,  the  trajectories  all  appear  to  converge  in  to  the  orbital  structure. 
The  controller  is  still  working  at  damping  large  perturbations  but  another  attractor 
has  appeared  in  the  phase  space.  Any  point  in  its  basin  of  attraction  will  result  in 
a  steady  state  at  that  non-linear  attractor.  To  gain  more  information  about  this 
attractor  will  require  some  other  non-linear  analysis  tools,  especially  the  Poincare 
map. 
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Fifure  4.23  Non-linear  ij.  to  ij_  for  200  orbit* 


4-2.4  Poincare  Map  To  observe  a  strange  attractor  with  a  Poincare  map  in 
a  space  higher  that  two  dimensions  requires  a  little  ingenuity.  Dr.  Moon  presented 
ideas  on  how  to  slice  up  a  three  dimensional  attractor  and  view  portions  of  it  at  a 
time  and  then  piece  together  the  whole  attractor.  The  real  difficulty  in  viewing  this 
attractor  is  the  fact  it  resides  in  a  four  dimensional  space.  This  can  be  achieved  in 
our  modal  coordinates  by  using  what  I  call  a  dual  torus,  Figure  4.25.  This  is  done 
by  stacking  two  doughnuts  on  top  of  one  another  and  running  the  time  coordinate 
as  the  angle  in  the  center  of  the  tori.  They  are  synchronized  to  the  same  time  axis. 
Then  r)i  vs  r/2  axis  is  at  the  center  of  the  first  torus  ring  and  r/:J  vs  ?/.,  axis  is  at  the 
center  of  the  second  torus  ring.  This  will  allow  us  to  display  all  four  dimensions  and 
take  a  time  slice  over  n  periods.  Each  system  state  is  represented  by  a  pair  of  points. 
Since  the  system  has  doubled  the  original  period,  start  the  map  with  a  4^  period. 

The  system  of  equations  are  integrated  for  many  orbits,  one  point  is  plotted 
in  each  coordinate  system  as  the  trajectories  pass  through  the  zero  angle  point  on 
the  tori.  Figures  4.26  and  4.27  show  the  results  of  2500  orbits.  There  are  four 
well  defined  regions  mapped  in  the  phase  space.  At  first  glance  they  do  not  appear 
to  be  fractal  nor  do  they  fit  any  of  the  other  cataloged  forms.  The  few  scattered 
points  off  of  the  shapes  are  the  initial  startup  disturbance,  before  it  settles  into  the 
attractor.  Let  us  take  a  series  of  slices  around  the  torus  and  examine  the  shape  over 
one  47T  period.  Ten  maps  are  shown  evenly  spaced  around  the  tori.  Figures  4.28, 
4.29,  4.30,  4.31,  4.32,  4.33,  4.34,  4.35,  4.36  and  4.37  show  the  attractor  in  the 
r]i  vs  t]2  space.  Figures  4.38,  4.39,  4.40,  4.41,  4.42,  4.43.  4.44,  4.45,  4.46  to 
4.47  show  the  r/3  vs  7?4  phase  space.  The  shape  of  the  attractor  seems  to  remain  in  a 
ring  type  shape  with  no  apparent  fractal  geometry.  One  thing  to  note,  examine  the 
shape  on  the  left  of  Figure  4.28  and  trace  it  through  the  4zr  cycle,  it  does  not  return 
back  to  the  same  region,  but  to  the  region  that  is  the  bottom  shape  in  Figure  4.28. 
Continuing  this  process  for  a  total  of  four  times  will  map  back  to  the  original  region 
on  the  left.  This  indicates  period  for  the  attractor  is  not  Air.  The  attractor  actually 
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has  a  period  of  167T.  Increasing  the  time  distance  around  the  tori  to  167T  and  making 
a  Poincare  map  of  the  data  can  be  seen  in  Figures  4.48  and  4.49.  Time  slices  could 
be  taken  to  examine  the  attractor  as  it  changes  in  time. 

4-3  Bifurcations  of  Cycles 

Where  did  this  attractor  come  from?  How  can  it  be  removed  or  compensation 
added?  These  questions  are  important  to  the  designer  and  operator  of  the  satellite. 
The  hint  lies  back  in  the  linearized  system.  The  period  of  the  system  has  doubled, 
this  is  usually  an  indicator  of  a  process  called  bifurcation.  A  process  called  a  Hopf 
bifurcation  has  been  examined  in  lower  order  systems  that  transitions  a  stable  node 
in  to  a  limit  cycle.  An  extension  of  this  idea  is  the  secondary  Hopf  bifurcation,  where 
a  limit  cycle  transitions  to  a  torus  [15:126-128].  In  our  case,  setting  the  gain  level 
beyond  the  bifurcation  point  has  caused  the  appearance  of  a  4-dimensional  torus 
around  the  equilibrium  point.  By  lowering  the  gain  below  the  bifurcation  point,  the 
torus  attractor  should  disappear. 

Using  the  starting  points  from  Table  4.5,  Figures  4.50  and  4.51  present 
typical  results  from  several  searches.  It  appears  that  the  equilibrium  point  is  now 
globally  stable  and  very  non-linear.  The  controller  is  effective  over  the  entire  phase 
space.  However  this  increased  stability  region  has  not  come  without  a  price.  The 
controller  now  takes  a  longer  time  to  dissipate  a  given  perturbation. 

By  slowly  increasing  the  gain  and  examining  the  phase  space,  the  appearance 
of  the  new  attractor  can  be  observed.  Figures  4.52,  4.53,  4.54,  4.55,  4.56,  4.57. 
4.58  and  4.59  show  the  formation  of  the  arractor.  This  provides  the  designer  with 
enough  information  to  specify  upper  gain  limits  for  the  controller.  With  the  gain  set 
at.  0.36  and  starting  at  fj  =  (2, 0. 0, 0} 7  and  examing  the  system  behavior  after  75 
orbits  we  find  the  limit  cycle  structure.  Figures  4.60  and  4.61. 

Also  a  two  dimensional  analogy  can  be  developed  that  provides  some  insight 
into  this  system,  figure  4.62.  This  system  behaves  similar  to  a  marble  in  a  bowl. 
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As  we  adjust  the  gain,  the  contour  of  the  bottom  of  the  bowl  changes.  With  a  gain 
setting  of  0.0,  the  system  exhibits  bounded  behavior  over  a  region  of  the  phase  space 
or  the  bottom  of  the  bowl  is  flat  and  has  no  preferred  equilibrium.  If  the  gain  is 
raised  to  0.2,  the  system  is  stable  over  the  entire  phase  space  and  is  attracted  to  only 
the  equilibrium  point  or  this  bowl  has  a  conventional  shape  with  trie  center  being  the 
lowest  point.  At  a  gain  of  0.35,  the  system  is  developing  a  limit  cycle  type  attractor 
or  the  bowl  now  has  a  near  flat  region  surrounding  the  center.  With  the  gain  at  0.4, 
the  phase  space  is  clearly  divided  into  two  regions  and  large  perturbations  will  not 
ret  urn  to  the  equilibrium  point  or  the  bowl  is  now  a  bowl  surrounded  by  a  ringshape. 

One  could  proceed  on  to  determine  exactly  what  type  of  attractor,  strange  or 
quasiperiodic,  is  appearing  with  the  bifurcation.  For  controller  design,  and  evaluation 
this  information  is  not  necessary  but  could  be  obtained  from  a  power  spectrum. 
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Table  4.5.  Starting  values  for  rji  axis  using  low  gain  controller 


m 

72 

73 

7-t 

0.5 

0.0 

iiiiiinrii 

1.0 

0.0 

0.0 

0.0 

1.5 

0.0 

0.0 

0.0 

2.0 

0.0 

0.0  I 

0.0  I 

3.0 

0.0 

ipjnii 

0.0 

5.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

n  e i hi  i  • 
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Figure  4.58  ij.  »«  rj-  for  100  orbits,  gain  =  0.35 


w; 


Figure  4.62.  Two  dimension  analogy  of  system 
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V.  Conclusions  and  Recommendations 


The  linearized  system  does  give  significant  insight  into  the  non-linear  problem. 
Myers  was  correct  in  picking  the  intersection  of  the  real  axis  as  being  the  most  stable 
setting  for  the  linearized  system  gain.  Since  this  is  where  the  bifurcation  occurs,  this 
implies  it  is  the  point  of  highest  gain  that  provides  global  stability.  However,  from 
the  non-linear  analysis,  we  can  recommend  gain  levels  up  to  0.36  and  still  have 
global  stability.  If  the  risk  of  large  perturbations  is  small,  even  highei  gains  could 
be  suggested  with  appropriate  limits.  The  non-linear  analysis  of  the  system  has  not 
produced  an  equation  that  will  tell  the  control  designer  what  value  of  gain  to  use. 
It  does  present  the  tools  required  to  determine  the  limits  of  a  controller  in  its  phase 
space  and  what  factors  can  effect  its  performance.  Each  controller  will  have  to  be 
evaluated  individually. 

Some  caution  needs  to  be  taken  since  quick  evaluations  around  the  linearized 
equilibrium  can  be  misleading.  Long  computer  runs  over  a  variety  of  starting  lo¬ 
cations  is  essential  in  evaluating  the  performance  of  a  control  system  design,  '[’lie 
uncontrolled  system  should  also  be  examined  closely  to  try  to  understand  any  non¬ 
linear  features.  Evaluation  using  the  full  non-linear  equations  of  motion  is  justified  in 
many  instances.  The  study  of  higher  order  systems  may  require  additional  computer 
resources. 

Sometimes  performance  must  be  traded  for  increased  range  of  stability.  This 
would  be  true  in  the  example  problem  where  one  would  have  to  set  the  gains  for 
longer  settling  times  to  be  prepared  for  the  case  the  satellite  should  be  bumped  very 
far  from  the  equilibrium.  One  could  not  risk  burning  up  the  controller  resources  and 
losing  a  satellite  just  to  provide  a  faster  correction. 

The  same  techniques  that  have  been  used  to  examine  this  scalar  controller 
could  be  applied  to  evaluate  other  controllers.  The  vector  controller  in  the  same 
satellite  is  an  obvious  next  step. 


Appendix  A.  Equations  of  Motion 


This  appendix  develops  the  equations  of  motion  of  a  spinning  inertially  sym¬ 
metric  satellite  in  an  elliptical  orbit  around  a  symmetrical  attracting  body.  The 
spin  axis  of  the  satellite  is  along  the  inertial  axis  of  symmetry  and  perpendicular  to 
the  orbit  plane.  The  equations  that  describe  the  attitude  of  the  satellite  are  pre¬ 
sented  and  they  are  linearized  about  an  equilibrium  point.  Both  the  linearized  and 
non-linear  attitude  equations  of  motion  are  developed. 

The  motion  and  stability  of  a  symmetrical  satellite  in  an  elliptical  orbit  has 
been  examined  previously.  Kane  and  Barba  developed  the  linearized  equations  that 
describe  the  attitude  motion  of  the  satellite.  This  section  is  taken  from  their  work 
[8]  and  the  expressions  for  the  non-linear  equations  of  motion  were  taken  from  the 
thesis  by  Captain  Dale  Shell  [14:2.23-2.25]. 


A.  I  Orbit  Equations 

Two  basic  equations  describe  the  trajectory  of  a  satellite  in  orbit  around  a 
spherically  symmetric  attracting  body.  These  equations  are  developed  in  any  intro¬ 
ductory  astrodynamics  text,  for  example  see  [16].  The  equations  are: 


r  —  n r  -f 


n2a3 


=  0 


(A.l) 


and, 

r2v  =  a2n(l  -  e2)1/2  (A.2) 

where  the  ”mean  motion”,  n,  is 

n  =  2*/T  (A. 3) 

and  T  is  the  period  of  the  orbit.  The  orbital  elements  r,  v,  a  and  e  are  defined  in 
f  igure  A.l.  r  and  v  are  the  polar  coordinates  of  the  center  of  mass  of  the  satellite. 


A-l 


A. 2  to  solve  for  v 


a2n( 1  —  e2) 
v  =  - o - 


2x1/2 


(A. 6) 


Substituting  this  into  equation  A.l  yields 


a4n2(l  —  e2)  n2a3 
r -  ' 


_3  r  2 

r  r 


H - =  0 


(A. 7, 


To  nondimensionalize  this  equation,  we  must  define  differentiation  with  respect  to 
the  nondimensional  time,  r.  All  derivatives  that  are  taken  with  respect  to  r  will  be 
shown  as  primes. 

dr  d  ( dr  \  d  . 

r  =  T,  =37UrJ  =SK<.)  =  (n« 

and 

r 


(A. 8) 


d  (  dr  \  d  ( dr  (  dr  \  \  H  2 

dt\dt)  dr  \dt  \dt ) ) 

'Sow  divide  equation  A. 7  by  n2a  to  get 


(A.9) 


n  a 


(A. 10) 


Then  substitute  in  equations  A. 8  and  A.9  this  reduces  to  the  nondimensional  form 


C"  +  72  + 


1  e2  -  1 


C  CJ 


=  0 


(A. 11) 


To  nondimensionalize  equation  A. 6  using  equations  A. 4  and  A. 5  note 

dv  du  dt  ( a2n{\  —  t2)xl2\  (  \ 


dr  dt  dr 


n 


this  reduces  to 


,  (1-e2)1'2 

v  =  - ^ - 

c2 


(A. 12) 


(A. 13) 
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and  differentiate  again  with  respect  to  r 


d  /(1-e2)1/2\ 

dr  V  C2  ) 


(A. 14) 


or 


// 


v 


2(l-e2)I/2., 

c3  c 


(A.  15) 


This  is  the  nondimensional  form  of  both  orbit  differential  equations.  The  initial 
conditions  are  arbitrarily  selected  such  that  the  satellite  is  at  the  perigee  when  t  = 
t  =  0.  It  then  follows  that 


r{ 0)  =  a(  1  -  e) 


(A. 16) 


and 


r(0)  =  0 


(A. 17) 


and  converting  it  to  the  nondimensional  form 


C(0)  =  1— e 


(A. 18) 


and 

C'(0)  =  0  (A.  19) 

It  is  now  possible  to  numerically  integrate  this  set  of  equations  to  find  the  path  of 
the  center  of  mass  of  the  satellite.  It  will  be  necessary  to  have  solutions  to  these  two 
equations  to  provide  inputs  to  the  satellite  attitude  equations. 

A. 2  Attitude  Equations 

The  orbit  equations  of  motion  assume  the  motion  of  the  center  of  mass  of  the 
satellite  is  unaffected  by  the  orientation  of  the  satellite.  The  attitude  equations  are 
not,  however,  independent  of  the  orbit  motion. 

Figure  A. 2  defines  an  orbital  reference  frame  A.  Ai  points  outward  along 
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Figure  A. 2.  Orbit  Reference  Frame 


the  orbit  radius,  A 2  is  perpendicular  to  A\  and  in  the  orbit  plane.  At  perigee  the 
satellite  is  moving  only  in  the  positive  A2  direction.  A3  is  perpendicular  to  .4!  and 
■42,  pointing  out  of  the  orbit  plane  to  form  a  right-hand  system. 

The  body  axis  frame  D  is  obtained  by  a  1-2-3  rotation  from  the  A  frame  as 
shown  in  Figure  A. 3.  The  three  angles,  &i,  02  and  $3  describe  the  magnitude  and 
direction  of  the  three  rotations.  <f>  is  the  angle  measured  from  the  .43  axis  to  the 
spin  axis  of  the  satellite,  D3.  Taking  advantage  of  the  symmetry  of  the  satellite,  the 
nodal  axis  C  can  be  used  to  simplify  the  analysis.  The  C  axis  is  obtained  through 
only  a  1-2  rotation.  This  leaves  us  one  rotation  away  from  the  body  axis,  but  the 
body  is  symmetric  along  this  axis.  To  develop  an  expression  for  the  angular  velocity 
of  the  satellite  in  the  nodal  coordinates  C,  we  start  by  writing  the  angular  velocity 
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of  the  reference  frame  A  with  respect  to  an  inertial  frame. 


UJA/I  =  ueM 


=  < 


0 

0 

v 


The  relative  rotation  rates  of  the  other  reference  frames  are 


t  \ 


0 1 

ojB/A  =  OxeB,  =  &xeAx  =  <j  0 

0 


uclB  =  02  es3  =  0i  ec2  =  s 


U>D/C  =  03ec3  =  03eo3  = 


0 

01  > 
0 

0 

0  } 

03 


(A. 20) 


(A.21) 


(A. 22) 


(A. 23) 


and  the  angular  velocity  of  the  satellite,  expressed  in  the  nodal  axes  is  the  sum 
of  these  components 

(A. 24) 
c 

To  add  these  components  of  the  angular  velocity,  they  must  be  expressed  in  the  same 
reference  frame  and  hence  we  transform  all  components  to  the  nodal  axis  C.  For 
more  information  on  nodal  axis  usage  can  be  found  in  [10]. 


jD'1  = 


u 


u>i 
id  2 

U>3 


K"}c  =  K'Tc  +  K"}fl  +  Lc*  K"}„  +  fco  K'T,  <*•»> 


A-T 


Where  the  transformation  matrices  are  defined  as 


Lcb  = 


Lb  a  = 


cos  $2  0  —sin  02 
0  1  0 

sin  $2  0  cos  02 

1  0  0 

0  cos  Q\  sin$! 

0  —sin  cos 0\ 


(A. 26) 


(A.  21 


cos  02  sin  0i  sin  02 

—  cos  0i  sin  02 

Lca  =  LcbLba  = 

0 

COS  0] 

sin  0! 

(A. 28) 

sin  0' 

>  —sin 0i  cos 02 

COS  0i  cos  02 

cos  03  sin  03 

0 

Lcd  = 

—  sin  03  cos  03 

0 

(A. 29) 

r 

o 

o 

1 

Then,  making  the  appropriate  substitutions  into  equation  A. 25,  all  components  of 
the  angular  velocity  are  transformed  into  the  nodal  reference  frame. 


{“D/'}c  = 


3  02  0  —  sin  02 


cos  02  sin  0i  sin  02 
0  cos  0i 


sin  02 

—  sin  0i 

cos  0 

’2 

COS  0 

cos  03 

sin  03 

0 

0 

-  sin  03 

cos  03 

0 

< 

0 

0 

0 

1 

03 

+ 


(A. 30) 
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After  performing  the  matrix  algebra,  this  reduces  to  the  final  expression  for  the 
angular  velocity  of  the  satellite  in  nodal  coordinates. 


K'%  = 


9\  cos  92  —  i>  cos  6i  sin  92 
62  +  v  sin  6 j 

9 1  sin  #2  +  v  cos  9\  cos  92  +  03 


(A. 31 ) 


By  defining  u?3  =  u>3  —  93  then  the  angular  velocity  can  be  expressed  as 


{-D/,}c = n + k,cl  - 


’  > 

U!X 

^1 

UJ2 

>  =  < 

UJ  2 

c 

t2’3  +  93 

.  4 

(A. 32) 


A  Iso, 


<i>3  =  9y  sin  92  +  v  cos  9i  cos  92 


(A. 33) 


This  expression  will  be  useful  in  the  development  of  Euler’s  moment  equations. 

The  angular  acceleration  can  be  found  by  differentiating  the  angular  velocity 
vector  with  respect  to  time  [9]. 


n = in 


(A. 34) 


where 


dt 


Sc 


’ 

'  > 

U>1 

Oi 

Cj  2 

’  =  ’ 

02 

W 3 

C 

a  3 

(A. 35) 
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and 


K"L*K"L  = 


'  ■> 

Wi 

u>i 

u?2 

>  X  ■ 

U)2 

.  =  < 

< 

c 

i  J 

c 

— U>iU>3  — 

U>iU,’2  —  U>iU>2 


•  ' 

0^u2 

►  =  < 

c 

0 

(A. 36) 

This  yields 


{“D,,}c  = 


—v  cos  Bi  sin  02  +  v  sin  sin  02  —  02  cos  Bl  cos  0 2)  + 
£  sin  #1  +  BB\  cos  B\  + 

^  u  cos  0X  cos  B2  —  0  [0 1  sin  B\  cos  &2  +  02  cos  Bx  sin  02)  + 

Bx  cos  82  —  6\B2  sin  B2 
02 


Bi  sin  02  +  0iO2  sin  02  +  0 3 


(A. 37) 


Having  solved  for  the  angular  acceleration,  use  the  angular  velocity  and  write  the 
angular  momentum  about  the  satellite  center  of  mass,  5. 


(Hs}C=U5)cK"}( 


(A. 38) 


Since  this  satellite  is  symmetric  about  the  spin  axis,  the  inertia  matrix  can  be  written 
as  a  matrix  of  constants 

A  0  0 

[Isle  =  0  A  0  (  (A. 39) 

0  0  c 

From  Euler’s  equation  of  motion,  the  change  in  angular  momentum  is  equal  to  the 
applied  moments. 

d-{Hs}c  =  {Ms}c  (A. 10) 


^{Hs}c+[uc"}{Hs}c  =  {Ms}c 


(A. 41) 
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where 


(A. 42) 


0  — UJ3  U>2 


—U>2  U)  2  0  J 

Since  the  inertia  matrix  is  constant,  equation  A. 41  can  be  written  as 
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(A. 43) 


This  reduces  to 


(  \  (  \ 


Aui\ 

0 

—U>3 

UJ2 

Au>\ 

Mi 

ACj  2 

•  + 

a>3 

0 

-UJX 

i 

Auj2 

>  =  < 

m2 

CC03 

— U>2 

U>2 

0 

Cu> 3 

m3 

si  j  v  s  \  ) 


Or  in  scalar  form 


(A. 44) 


Mi  —  Aujj  —  Auj 2^3  A  Cu>2u?3 

(A. 45) 

M2  —  AcL>2  -b  Alu\LJ3  —  Cu) JW3 

(A. 46) 

M3  —  Cu3  —  AuJiU)2  T  AuJ] ui2 

(A. 47) 

Replacing  u3  with  u3  —  03  this  becomes 

M\  —  A 0J\  —  A0J2OJ3  “f"  AuJ 2^3  “1“  CuJ2^3 

(A. 48) 

A/2  — 3  ACj2  "("  Aujiljs  —  AuJxO^  —  Cu}\u?$ 

(A. 49) 

•s 

0 

II 

(A. 50) 
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Dividing  the  first  two  moment  equations  by  A,  and  the  last  by  C . 


Mx 

A 

C  —  A  • 

(A. 51) 

—  u.'l  + 

Ul2  UJ3  +  UJ2V3 

A 

m2 

C  -  A  • 

(A. 52) 

A 

=  Ul2  — 

UJ\UJ3  U)\t)2 

A 

m3 

C 

=  <I?3 

(A. 53) 

The  moments  from  gravitational  potentials  are  given  in  [8].  These  values  are 


Mi  =  0  (A.  54) 

A/2  =  —  3n2(~3  (C  —  .4)  sin  &2  cos  02  (A. 55) 

A/3  =  0  (A. 56) 


Introducing  the  inertial  ratio  I\  that  is  defined  as 

the  three  moment  equations  can  be  written  as 


UJ\  "1"  A  U?2U^3  4"  lc?2$3  =  0 

(A. 58) 

u>2  —  Kuj\u>3  —  —  —  3n2(”~3A'sin^2COS^2 

(A. 59) 

0 

II 

n 

•3 

(A. 60) 

QTi  K U>2W3  +  U/203  —  0 

(A. 61) 

a2  —  Kuxu>3  —  uq#3  =  —  3n2£-3 K  sin  02  cos  02 

(A. 62) 

q3  =  0 

(A. 63) 

where  {a}  is  given  in  equation  A. 37  and  {u>}  is  from  equation  A. 31. 
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Now  we  consider  motion  of  the  satellite  about  the  equilibrium  point  0,  =  02  = 
0.  The  value  of  03  can  be  solved  for  in  terms  of  u>3.  Using  small  angle  conditions 
around  the  equilibrium  point 


r 

U>1 

'  > 

0i  —  i/02 

>  =  i 

0 2  +  1/0 1 

OJ3 

c 

1/  +  03 

>  j 

Integrating  the  equation  for  u;3 


—  u;3f  =  $3  +  03(0)  +  1/  +  t/(0) 


(A.  64) 


(A. 65) 


With  the  initial  conditions  defined  as  03(O)  =  i/(0)  =  0  solve  for  03. 


03  =  u/3 1  —  v  (A. 66) 

where  u;3  =  n(l  +  a),  a  is  the  spin  rate  of  the  satellite,  |u/D/c'|,  and  is  a  constant 
value,  n  is  defined  in  equation  A. 3.  For  these  conditions,  the  accelerations  and 
moments  are  zero.  I  herefore  this  represents  a  valid  equilibrium  point,  with  the  spin 
axis  perpendicular  to  the  orbit  plane.  .  The  next  step  is  to  transform  our  linearized 
differential  equations  of  motion  into  a  form  where  the  stability  of  this  equilibrium 
point  can  be  evaluated  using  Floquet  Theory.  The  desired  form  is 

x'  =  Ax  (A. 67) 
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where  x  represents  the  vector  of  state  variables. 


{*}  = 


Oi 

02 

01 

02 


(A. 68) 


The  first  two  rows  of  the  A  matrix  are  fairly  simple  with  a  one  in  the  appropriate 
column  and  the  rest  of  the  elements  zero.  The  last  two  rows  can  be  found  by 
differentiating  equation  A. 64. 


(It  t  Jc 


' 

’  ' 

u;i 

U>2 

►  -  < 

a2 

>  =  < 

U)3 

.  > 

C 

C*3 

V.  J 

C 

Q\  —  r/02  —  0O2 
&2  T  v9\  T  00\ 

0  0  3 


(A. 69) 


Substituting  these  values  for  a  and  the  values  for  u>  from  equation  A. 64  into  equa¬ 
tions  A. 62  and  A. 63  the  result  is  two  linearized  equations  for  the  motion  in  the 
region  near  the  equilibrium  point. 

0i  —  002  —  00 2  +  A  ^ 02  +  09\^j  (03  +  T  (&2  T  00\j  O3  =  0  (A.  (0) 

also  applying  the  small  angle  assumptions  to  the  second  equation 

02  +  00i  +  09\  —  K  ^0\  —  002^j  ( 03  +  0^j  —  (Oi  —  002 'j  03  =  —3n2  A  02^  1  (A. 71) 

To  put  these  in  the  proper  format  we  must  apply  the  chain  rule  and  convert  them  to 
a  form  with  the  differentiation  is  with  respect  to  r.  Then  solve  for  0"  and  02.  This 
results  in 

0"  =  -i/'  (0'  +  AV  +  K0'3)  0,  +  i/"02  +  («/'  -  0'  -  K  (!/'  +  0' ))  0'  (A. 72) 
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02  =  -»'%  +  [-«/'  {0'3  +  AV  +  I\03)  -  3A'C3]  O2-W  -O3-  A  w  +  0'3))  0\  (A. 73) 


and  u"  can  be  replaced  by  the  values  from  equations  A. 13  and  A. 14.  0'3  can  be 
replaced  by  noting  that  a3  =  0,  then  from  equation  A. 64 

v  +  03  =  constant  (A. 74) 

and  letting  the  constant  be  a  +  1,  then  0'3  can  be  expressed  as 

0'3  =  (a  +  l)-r2(l-e2)l/2  (A. 75) 

Now  the  equations  can  be  expressed  in  matrix  form  in  terms  satellite  inertia  ratio 
(A')  ,  satellite  spin  rate  (a)  and  orbital  eccentricity  (C  and  (').  If  the  spin  rate  and 
inertia  ratio  of  the  satellite  are  held  fixed,  increased  eccentricity  will  tend  to  desta¬ 
bilize  the  satellite  attitude. 


where: 


A(t)  = 


0  0 

0  0 

fl31  <*32 
a4j  a4 2 


1  0 

0  1 

0  a34 

a43  0 


a31  =  (1  -  e2)C"4  -  (1  +  a)(l  +  I<)(1  -  e2),/2C~2 
032  =  -041  =  -2(1  -  e2)l/2C'C"3 

034  =  -043  =  2(1  -  -  (1  +  a)(l  +  K) 

042  =  a31  —  3A'C-3 


(A.76) 


It  should  be  noted  all  the  elements  of  A  are  periodic  with  a  period  of  one  orbit,  T. 
This  is  the  proper  form  needed  to  apply  Floquet  Theory.  Now  the  stability  can  be 
determined  in  the  region  of  an  equilibrium  point. 
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A. 3  Non-linear  Equations  of  Motion 

All  of  the  previous  development  has  been  based  on  the  linearized  equations  of 
motion  around  an  equilibrium  point.  To  fully  evaluate  the  system  requires  its  non¬ 
linear  behavior  to  be  examined  in  the  phase  space.  To  do  this  requires  the  non-linear 
equations  of  motion.  The  expressions  for  these  can  be  found  by  replacing  u?i,  u?2,  on 
and  a 2  in  equations  A. 62  and  A. 63  with  the  expressions  from  equations  A. 31  and 
A. 37.  is  still  a  constant,  a  4-  1,  that  is  the  spin  rate  about  the  axis  of  symmetry. 
Using  this  value  and  solving  equation  A. 31,  to  remove  03  from  the  equations  yields. 

^  =  a  +  1  —  Q\  sin  &2  —  v  cos  0t  cos  d2  (A. 77) 

Making  the  appropriate  substitutions  into  equation  A. 62  yields. 

O  =  0X  cos  02  —  0i02  sin  02  —  v  cos  0!  sin  02  -f  i/0\  sin  0!  sin  02  —  i>02  cos  0 1  cos  02  + 

I\  (02  +  I>sin0i)(a  +  1)  +  (02  +  i>sin0i)03  (ATS) 

Converting  the  derivative  from  clt  to  dr,  replacing  03  then  solving  for  0" 

0"  =  29[92  tan  02  +  2 v'9'2  cos  0[  —  ( K  +  l)(a  +  1  )(9'2  +  v  sin  0i )  sec  02 

+  o"  cos  0i  tan  02  +  v'  v  sin  0i  cos  0t  (A. 79) 

And  in  a  similar  manner  solve  for  92 

02  =  (K  +  1)(q  +  l)(0j  cos  02  —  v  cos  0i  sin02)  —  v"  sin0i  —  '2v9\  cos0i  cos2  02  — 
(0j2  —  vn  cos2  0i  +  3A'C-3)  s'n  02  cos  02  (A. 80) 

Expressions  for  i/  and  v"  are  in  equations  A. 6  and  A. 7.  Numerical  integration  of 
this  equation  set  will  allow  us  to  determine  the  motion  of  the  satellite  described  by 
the  non-linear  equations  of  motion. 
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